A synthesis of dryland restoration techniques.

Purpose

To quantitatively examine the efficacy of vegetation restoration in drylands globally.

Questions

1. What were the primary restoration goals as reported by primary authors?
   
2. How much variation was there in the techniques tested and how long were experiments monitored and tested?
   
3. How effective are dryland restoration techniques for vegetation?

Step 1. Search

#general library loads
library(tidyverse)

#functions
se <- function(x){
  sd(x)/sqrt(length(x))
}

#study data####
studies <- read_csv("data/studies.csv")

#studies included
evidence <- studies %>%
  filter(exclude == "no")

dist<- ggplot(evidence, aes(disturbance, fill = paradigm)) +
  geom_bar(na.rm = TRUE) +
  coord_flip() +
  scale_fill_brewer(palette = "Paired") +
  labs(y = "frequency")+
  theme(panel.background = element_rect(fill="white", colour="gray50",linetype="solid"))+
  theme(panel.grid.minor = element_line(colour="gray80", linetype="dashed"))

#ggsave(file="dist.svg", width=9, height=5)

ggplot(evidence, aes(intervention, fill = paradigm)) +
  geom_bar(na.rm = TRUE) +
  coord_flip() +
  scale_fill_brewer(palette = "Paired") +
  labs(y = "frequency")

#paradigm
derived.evidence <- evidence %>%
  group_by(technique, data, region, disturbance, goal, paradigm) %>% summarise(n = n())

#active-passive split
m <- glm(n~paradigm, family = poisson, derived.evidence)
anova(m, test="Chisq")

## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: n
## 
## Terms added sequentially (first to last)
## 
## 
##          Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL                       167     9.9147         
## paradigm  1 0.045115       166     9.8696   0.8318

#region
m1 <- glm(n~paradigm*region, family = poisson, derived.evidence)
#m1
#summary(m1)
anova(m1, test="Chisq")

## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: n
## 
## Terms added sequentially (first to last)
## 
## 
##                 Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL                              167     9.9147         
## paradigm         1 0.045115       166     9.8696   0.8318
## region           6 0.301367       160     9.5682   0.9995
## paradigm:region  6 0.213627       154     9.3546   0.9998

#outcome
m2 <- glm(n~paradigm*goal, family = poisson, derived.evidence)
#m1
#summary(m1)
anova(m2, test="Chisq")

## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: n
## 
## Terms added sequentially (first to last)
## 
## 
##               Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL                            167     9.9147         
## paradigm       1 0.045115       166     9.8696   0.8318
## goal           6 0.240941       160     9.6287   0.9997
## paradigm:goal  4 0.301480       156     9.3272   0.9897

Step 2. Sort

library(PRISMAstatement)
prisma(found = 1504,
       found_other = 5,
       no_dupes = 1039, 
       screened = 1039, 
       screen_exclusions = 861, 
       full_text = 178,
       full_text_exclusions = 101, 
       qualitative = 77, 
       quantitative = 40,
       width = 800, height = 800)

#pdf(prisma, "prisma.pdf")

Step 3. Synthesize

#all data includes non-relevant and inc search term studies
data_all <- read_csv("data/data_all.csv")

#data from ag & grazing studies that examined restoration in drylands
data <- data_all %>% 
  filter(disturbance %in% c("agriculture","grazing")) %>% 
  filter(!notes %in% "couldnt extract data") %>%
  mutate(lrr = log(mean.t/mean.c), var.es = ((sd.t^2/(n.t*mean.t^2)) + (sd.c^2/(n.c*mean.c^2)))) %>%
  filter(!is.na(lrr)) %>%
  filter(!is.na(var.es)) %>%
  filter(!is.na(n.t)) %>%
  filter(!is.na(p)) %>%
  filter(!is.na(intervention)) %>%
  filter(is.finite(lrr)) %>%
  filter(!is.na(exp.length)) %>%
  filter(!is.na(MAP)) %>%
  filter(!is.na(aridity.index))

#write cleaned data for provenace and more rapid reuse
#write_csv(data, "data/data.csv")


#totable <- data %>% group_by(paradigm, intervention, technique, outcome) %>% count()
#write.csv(totable,"totable.csv")

Step 4. Summarize

#evidence map####
require(maps)
world<-map_data("world")
map<-ggplot() + geom_polygon(data=world, fill="gray80", aes(x=long, y=lat, group=group))

#map + geom_point(data=data, aes(x=long, y=lat, color = paradigm), size=2)+ 
  #theme(panel.background = element_rect(fill="white", colour="gray50",linetype="solid"))+
  #theme(panel.grid.minor = element_line(colour="gray80", #linetype="dashed"))+
  #scale_color_brewer(palette = "Paired") +
  #theme(legend.position = c(0.1, 0.3)) +
  #guides(fill=guide_legend(title=NULL)) +
  #theme(legend.background = element_rect(fill="white", 
                                  #size=0.5, linetype="solid", #colour="gray70"))+
  #labs(x = "longitude", y = "latitude", color = "")


map + geom_point(data=data_all, aes(x=long, y=lat, color = paradigm), size=2)+ 
  theme(panel.background = element_rect(fill="white", colour="gray50",linetype="solid"))+
  #theme(panel.grid.minor = element_line(colour="gray80", linetype="dashed"))+
  scale_color_brewer(palette = "Set1") +
  theme(legend.position = c(0.1, 0.3)) +
  guides(fill=guide_legend(title=NULL)) +
  theme(legend.background = element_rect(fill="white", 
                                  size=0.5, linetype="solid", colour="gray70"))+
  labs(x = "longitude", y = "latitude", color = "")

#ggsave(file="map.svg", width=12.7, height=8)

Step 5. Models

#meta####
library(meta)
#active-passive differences####
m1 <- metagen(lrr, var.es, studlab = ID, comb.fixed = FALSE, byvar = paradigm, data = data)
summary(m1)

## Number of studies combined: k = 1460
## 
##                                       95%-CI     z  p-value
## Random effects model 0.0766 [0.0654; 0.0879] 13.38 < 0.0001
## 
## Quantifying heterogeneity:
## tau^2 = 0.0450; H = 184527270028.74 [184527270027.80; 184527270029.68]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Quantifying residual heterogeneity:
## H = 184583923516.49 [184583923515.55; 184583923517.44]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Test of heterogeneity:
##                              Q d.f. p-value
##  49679407227636228416848646.00 1459       0
## 
## Results for subgroups (random effects model):
##                       k                     95%-CI
## paradigm = active  1102  0.2184 [ 0.2055;  0.2314]
## paradigm = passive  358 -0.3413 [-0.3753; -0.3073]
##                                                Q  tau^2    I^2
## paradigm = active  49671693556322550580400868.00 0.0450 100.0%
## paradigm = passive     4152232320954931871684.00 0.1047 100.0%
## 
## Test for subgroup differences (random effects model):
##                       Q d.f.  p-value
## Between groups   911.23    1 < 0.0001
## 
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2

#funnel(m1)
#radial(m1)
#forest(m1, layout = "JAMA", study.results = FALSE)

#t-tests if different from 0
tmu <- function(x){t.test(x, mu = 0, paired = FALSE, var.equal=FALSE, conf.level = 0.95)
  }

data %>%
  split(.$paradigm) %>%
  purrr::map(~tmu(.$lrr)) #note this uses arithmetic means not estimated means from random effect models

## $active
## 
##  One Sample t-test
## 
## data:  x
## t = 7.6083, df = 1101, p-value = 5.943e-14
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.2069479 0.3507825
## sample estimates:
## mean of x 
## 0.2788652 
## 
## 
## $passive
## 
##  One Sample t-test
## 
## data:  x
## t = -7.5438, df = 357, p-value = 3.824e-13
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.4408271 -0.2585133
## sample estimates:
##  mean of x 
## -0.3496702

#metareg(m1, ~aridity.index+exp.length) #covariates and additive
mr1 <- metareg(m1, ~aridity.index*exp.length)#interaction term
mr1

## 
## Mixed-Effects Model (k = 1460; tau^2 estimator: DL)
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0440 (SE = 0.0167)
## tau (square root of estimated tau^2 value):             0.2098
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability):   33306623692228221992660.00
## R^2 (amount of heterogeneity accounted for):            2.22%
## 
## Test for Residual Heterogeneity:
## QE(df = 1456) = 48494444095884294040600408.0000, p-val < .0001
## 
## Test of Moderators (coefficients 2:4):
## QM(df = 3) = 5927.5166, p-val < .0001
## 
## Model Results:
## 
##                           estimate      se      zval    pval    ci.lb 
## intrcpt                     0.4843  0.0171   28.3215  <.0001   0.4508 
## aridity.index              -0.0142  0.0008  -18.3338  <.0001  -0.0158 
## exp.length                  0.0025  0.0001   35.5749  <.0001   0.0024 
## aridity.index:exp.length   -0.0002  0.0000  -48.9045  <.0001  -0.0002 
##                             ci.ub 
## intrcpt                    0.5178  *** 
## aridity.index             -0.0127  *** 
## exp.length                 0.0027  *** 
## aridity.index:exp.length  -0.0001  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(mr1)

#interventions####
#active
m2 <- metagen(lrr, var.es, studlab = ID, byvar = intervention, comb.fixed=FALSE, subset = paradigm == "active", data = data)
summary(m2)

## Number of studies combined: k = 1102
## 
##                                       95%-CI     z  p-value
## Random effects model 0.2184 [0.2055; 0.2314] 33.02 < 0.0001
## 
## Quantifying heterogeneity:
## tau^2 = 0.0450; H = 212403086959.62 [212403086958.53; 212403086960.71]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Quantifying residual heterogeneity:
## H = 210054235083.87 [210054235082.78; 210054235084.96]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Test of heterogeneity:
##                              Q d.f. p-value
##  49671693556322550580400868.00 1101       0
## 
## Results for subgroups (random effects model):
##                                 k                  95%-CI
## intervention = vegetation     779 0.1845 [0.1694; 0.1996]
## intervention = soil           248 0.3128 [0.2990; 0.3265]
## intervention = water addition  75 0.6409 [0.5539; 0.7279]
##                                                           Q  tau^2    I^2
## intervention = vegetation     48393423300974332080088226.00 0.0443 100.0%
## intervention = soil              97513761686148203152860.00 0.0111 100.0%
## intervention = water addition                      31132.18 0.1047  99.8%
## 
## Test for subgroup differences (random effects model):
##                       Q d.f.  p-value
## Between groups   226.58    2 < 0.0001
## 
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2

#funnel(m2)
#radial(m2)
#metabias(m2)
#forest(m2, layout = "JAMA", study.results = FALSE)

#metareg(m2, ~ aridity.index + exp.length)
mr2 <- metareg(m2, ~ aridity.index*exp.length)
summary(mr2)

## 
## Mixed-Effects Model (k = 1102; tau^2 estimator: DL)
## 
##     logLik    deviance         AIC         BIC        AICc 
## -8260.1776  38647.4255  16530.3553  16555.3797  16530.4100   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0440 (SE = 0.0167)
## tau (square root of estimated tau^2 value):             0.2098
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability):   44155839198290107172246.00
## R^2 (amount of heterogeneity accounted for):            2.23%
## 
## Test for Residual Heterogeneity:
## QE(df = 1098) = 48483111439722536500202208.0000, p-val < .0001
## 
## Test of Moderators (coefficients 2:4):
## QM(df = 3) = 2277.4091, p-val < .0001
## 
## Model Results:
## 
##                           estimate      se      zval    pval    ci.lb 
## intrcpt                     0.5600  0.0193   29.0873  <.0001   0.5222 
## aridity.index              -0.0145  0.0009  -15.8801  <.0001  -0.0163 
## exp.length                  0.0033  0.0001   26.6964  <.0001   0.0030 
## aridity.index:exp.length   -0.0002  0.0000  -20.9978  <.0001  -0.0003 
##                             ci.ub 
## intrcpt                    0.5977  *** 
## aridity.index             -0.0127  *** 
## exp.length                 0.0035  *** 
## aridity.index:exp.length  -0.0002  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(mr2)

mr2.1 <- metareg(m2, ~ aridity.index)
bubble(mr2.1)

mr2.2 <- metareg(m2, ~ exp.length)
bubble(mr2.2)

#passive
m3 <- metagen(lrr, var.es, studlab = ID, byvar = intervention, subset = paradigm == "passive", comb.fixed=FALSE, data = data)
summary(m3)

## Number of studies combined: k = 358
## 
##                                          95%-CI      z  p-value
## Random effects model -0.3413 [-0.3753; -0.3073] -19.70 < 0.0001
## 
## Quantifying heterogeneity:
## tau^2 = 0.1047; H = 3410410951.75 [3410410950.14; 3410410953.36]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Quantifying residual heterogeneity:
## H = 3420000814.26 [3420000812.64; 3420000815.87]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Test of heterogeneity:
##                          Q d.f. p-value
##  4152232320954931871684.00  357       0
## 
## Results for subgroups (random effects model):
##                                    k                     95%-CI
## intervention = vegetation        125  0.2654 [ 0.2067;  0.3241]
## intervention = grazing exclusion  29  0.1351 [ 0.0270;  0.2431]
## intervention = soil              204 -0.7583 [-0.8196; -0.6970]
##                                                          Q  tau^2    I^2
## intervention = vegetation        4152209524073903423668.00 0.1047 100.0%
## intervention = grazing exclusion              238316232.18 0.0881 100.0%
## intervention = soil                   14453104123321616.00 0.1990 100.0%
## 
## Test for subgroup differences (random effects model):
##                       Q d.f.  p-value
## Between groups   595.91    2 < 0.0001
## 
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2

#funnel(m3)
#radial(m3)
#metabias(m3, method = "linreg")

#metareg(m3, ~ aridity.index + exp.length)
mr3 <- metareg(m3, ~ aridity.index*exp.length)
summary(mr3) 

## 
## Mixed-Effects Model (k = 358; tau^2 estimator: DL)
## 
##    logLik   deviance        AIC        BIC       AICc 
## -675.8601  5439.3040  1361.7201  1381.1228  1361.8906   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.1047 (SE = 0.0876)
## tau (square root of estimated tau^2 value):             0.3236
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability):   11729120705162215424.00
## R^2 (amount of heterogeneity accounted for):            0.00%
## 
## Test for Residual Heterogeneity:
## QE(df = 354) = 4152108729627424325662.0000, p-val < .0001
## 
## Test of Moderators (coefficients 2:4):
## QM(df = 3) = 882.1252, p-val < .0001
## 
## Model Results:
## 
##                           estimate      se     zval    pval    ci.lb 
## intrcpt                     0.1386  0.0636   2.1770  0.0295   0.0138 
## aridity.index               0.0042  0.0029   1.4244  0.1543  -0.0016 
## exp.length                  0.0106  0.0015   6.8767  <.0001   0.0076 
## aridity.index:exp.length   -0.0005  0.0001  -8.1677  <.0001  -0.0006 
##                             ci.ub 
## intrcpt                    0.2633    * 
## aridity.index              0.0099      
## exp.length                 0.0136  *** 
## aridity.index:exp.length  -0.0004  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(mr3)

mr3.1 <- metareg(m3, ~ aridity.index)
bubble(mr3.1)

mr3.2 <- metareg(m3, ~ exp.length)
bubble(mr3.2)

#outcomes
#active
m4 <- metagen(lrr, var.es, studlab = ID, byvar = outcome, subset = paradigm == "active", comb.fixed=FALSE, data = data)
summary(m4)

## Number of studies combined: k = 1102
## 
##                                       95%-CI     z  p-value
## Random effects model 0.2184 [0.2055; 0.2314] 33.02 < 0.0001
## 
## Quantifying heterogeneity:
## tau^2 = 0.0450; H = 212403086959.62 [212403086958.53; 212403086960.71]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Quantifying residual heterogeneity:
## H = 210149866430.08 [210149866428.99; 210149866431.17]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Test of heterogeneity:
##                              Q d.f. p-value
##  49671693556322550580400868.00 1101       0
## 
## Results for subgroups (random effects model):
##                     k                     95%-CI
## outcome = soil    249  0.2204 [ 0.1558;  0.2849]
## outcome = plants  305  0.5071 [ 0.4936;  0.5206]
## outcome = animals  24 -0.1152 [-0.1155; -0.1148]
## outcome = habitat 524  0.0621 [ 0.0437;  0.0804]
##                                               Q   tau^2    I^2
## outcome = soil                35077220764051.67  0.2656 100.0%
## outcome = plants     97513760543782884872406.00  0.0111 100.0%
## outcome = animals                     541696.64 <0.0001 100.0%
## outcome = habitat 48393423303339003636244404.00  0.0443 100.0%
## 
## Test for subgroup differences (random effects model):
##                        Q d.f. p-value
## Between groups   8647.81    3       0
## 
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2

#metabias(m)
mr4 <- metareg(m4, ~aridity.index*exp.length)
mr4

## 
## Mixed-Effects Model (k = 1102; tau^2 estimator: DL)
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0440 (SE = 0.0167)
## tau (square root of estimated tau^2 value):             0.2098
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability):   44155839198290107172246.00
## R^2 (amount of heterogeneity accounted for):            2.23%
## 
## Test for Residual Heterogeneity:
## QE(df = 1098) = 48483111439722536500202208.0000, p-val < .0001
## 
## Test of Moderators (coefficients 2:4):
## QM(df = 3) = 2277.4091, p-val < .0001
## 
## Model Results:
## 
##                           estimate      se      zval    pval    ci.lb 
## intrcpt                     0.5600  0.0193   29.0873  <.0001   0.5222 
## aridity.index              -0.0145  0.0009  -15.8801  <.0001  -0.0163 
## exp.length                  0.0033  0.0001   26.6964  <.0001   0.0030 
## aridity.index:exp.length   -0.0002  0.0000  -20.9978  <.0001  -0.0003 
##                             ci.ub 
## intrcpt                    0.5977  *** 
## aridity.index             -0.0127  *** 
## exp.length                 0.0035  *** 
## aridity.index:exp.length  -0.0002  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(mr4)

#passive
m5 <- metagen(lrr, var.es, studlab = ID, byvar = outcome, subset = paradigm == "passive", comb.fixed=FALSE, data = data)
summary(m5)

## Number of studies combined: k = 358
## 
##                                          95%-CI      z  p-value
## Random effects model -0.3413 [-0.3753; -0.3073] -19.70 < 0.0001
## 
## Quantifying heterogeneity:
## tau^2 = 0.1047; H = 3410410951.75 [3410410950.14; 3410410953.36]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Quantifying residual heterogeneity:
## H = 3419999089.05 [3419999087.44; 3419999090.67]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Test of heterogeneity:
##                          Q d.f. p-value
##  4152232320954931871684.00  357       0
## 
## Results for subgroups (random effects model):
##                     k                     95%-CI                         Q
## outcome = habitat 104  0.1605 [ 0.0964;  0.2246] 4152172019980636258444.00
## outcome = plants   50  0.4438 [ 0.0345;  0.8532]      33314950066906688.00
## outcome = soil    204 -0.7583 [-0.8196; -0.6970]      14453104123321616.00
##                    tau^2    I^2
## outcome = habitat 0.1047 100.0%
## outcome = plants  2.1620 100.0%
## outcome = soil    0.1990 100.0%
## 
## Test for subgroup differences (random effects model):
##                       Q d.f.  p-value
## Between groups   425.72    2 < 0.0001
## 
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2

#metabias(m)
mr5 <- metareg(m5, ~ aridity.index*exp.length)
mr5

## 
## Mixed-Effects Model (k = 358; tau^2 estimator: DL)
## 
## tau^2 (estimated amount of residual heterogeneity):     0.1047 (SE = 0.0876)
## tau (square root of estimated tau^2 value):             0.3236
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability):   11729120705162215424.00
## R^2 (amount of heterogeneity accounted for):            0.00%
## 
## Test for Residual Heterogeneity:
## QE(df = 354) = 4152108729627424325662.0000, p-val < .0001
## 
## Test of Moderators (coefficients 2:4):
## QM(df = 3) = 882.1252, p-val < .0001
## 
## Model Results:
## 
##                           estimate      se     zval    pval    ci.lb 
## intrcpt                     0.1386  0.0636   2.1770  0.0295   0.0138 
## aridity.index               0.0042  0.0029   1.4244  0.1543  -0.0016 
## exp.length                  0.0106  0.0015   6.8767  <.0001   0.0076 
## aridity.index:exp.length   -0.0005  0.0001  -8.1677  <.0001  -0.0006 
##                             ci.ub 
## intrcpt                    0.2633    * 
## aridity.index              0.0099      
## exp.length                 0.0136  *** 
## aridity.index:exp.length  -0.0004  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(mr5)

detach(package:meta, unload = TRUE)

#library(metafor)
#must use dev version of metafor (https://wviechtb.github.io/metafor/#installation)

#data<-escalc(measure="ROM",m1i=mean.t,m2i=mean.c,sd1i=sd.t,sd2i=sd.c,n1i=n.t,n2i=n.c, #data=mydata,var.names=c("LRR","LRR_var"),digits=4)

#data <- data %>%
#  filter(!is.na(LRR)) %>%
#  filter(!is.na(LRR_var)) %>%
#  filter(!is.na(n.t)) %>%
#  filter(!is.na(p)) %>%
#  filter(!is.na(intervention)) %>%
#  filter(is.finite(lrr)) %>%
#  filter(!is.na(exp.length)) %>% 
#  filter(!is.na(aridity.index))

#mod.1 <- rma(yi=lrr, vi=var.es, mods = ~paradigm, data = data)
#summary(mod.1)
#forest(mod.1)

#interventions
#mod.2 <- rma(lrr, var.es, slab= ID, mods=  ~intervention+aridity.index*exp.length -1, data = data, subset = paradigm == "active")
#summary(mod.2)
#forest(mod.3, slab= "study.ID")

#mod.3 <- rma(lrr, var.es, slab= ID, mods=  ~intervention+aridity.index*exp.length -1, data = data, subset = paradigm == "passive")
#summary(mod.3)

#outcomes
#mod.4 <- rma(lrr, var.es, slab= ID, mods=  ~outcome+aridity.index*exp.length -1, data = data, subset = paradigm == "active")
#summary(mod.4)

#mod.5 <- rma(lrr, var.es, slab= ID, mods=  ~outcome+aridity.index*exp.length -1, data = data, subset = paradigm == "passive")
#summary(mod.5)

#detach(package:metafor, unload = TRUE)

Viz

#ggplot uses geometric means and random means reported in models are harmonic

ggplot(data, aes(paradigm, lrr, color = intervention)) +
  ylim(c(-2,2)) +
  geom_boxplot() + 
  labs(x = "", y = "lrr", color = "") +
  coord_flip() +
  geom_hline(yintercept = 0, colour="grey", linetype = "longdash")+
  theme(axis.text.x=element_text(face="bold"),
        axis.text.y=element_text(face="bold"),
        axis.title=element_text(size=12,face="bold"),
        strip.text.y = element_text(hjust=0,vjust = 1,angle=180,face="bold")) +
  scale_color_brewer(palette = "Set1")

ggplot(data, aes(paradigm, lrr, fill = intervention)) +
  ylim(c(-2,2)) +
  geom_violin() + 
  labs(x = "", y = "lrr", fill = "") +
  coord_flip() +
  geom_hline(yintercept = 0, colour="grey", linetype = "longdash")+
  theme(axis.text.x=element_text(face="bold"),
        axis.text.y=element_text(face="bold"),
        axis.title=element_text(size=12,face="bold"),
        strip.text.y = element_text(hjust=0,vjust = 1,angle=180,face="bold")) +
  scale_fill_brewer(palette = "Set1")

ggplot(data, aes(aridity.index, lrr, color = intervention)) + 
  geom_point(position = position_dodge(width = 0.5)) +
  facet_wrap(~paradigm) +
  scale_color_brewer(palette = "Set1") + 
  labs(x = "aridity", y = "lrr", color = "")

ggplot(data, aes(exp.length, lrr, color = intervention)) + 
  geom_point(position = position_dodge(width = 0.5)) +
  facet_wrap(~paradigm) +
  scale_color_brewer(palette = "Set1") + 
  labs(x = "length of experiment", y = "lrr", color = "")

#random model outputs
models <- read_csv("data/meta_outputs.csv")
ggplot(models, aes(paradigm, lrr, color = intervention)) +
  ylim(c(-2,2)) +
  geom_point(position = position_dodge(width = 0.5), size=2) + 
  labs(x = "", y = "log response ratio", color = "") +
  coord_flip() +
  geom_errorbar(aes(ymin=lower, ymax=upper), size= 0.7, width=0.2, position = position_dodge(width = 0.5)) +
  geom_hline(yintercept = 0, colour="grey", linetype = "longdash")+
  theme(axis.text.x=element_text(face="bold"),
        axis.text.y=element_text(face="bold"),
        axis.title=element_text(size=12,face="bold"),
        strip.text.y = element_text(hjust=0,vjust = 1,angle=180,face="bold")) +
  theme_bw()+
  scale_color_brewer(palette = "Set1")

ggplot(models, aes(paradigm, lrr, color = intervention))+
  ylim(c(-1,1)) +
  geom_point(position = position_dodge(width = 0.5), size=2) + 
scale_x_discrete(limits=c("passive","active"))+
  labs(x = "", y = "log response ratio", color = "") +
  coord_flip() +
  geom_text(data=models, mapping=aes(x=paradigm, y=lrr, label=k, group=intervention),  position = position_dodge(width = 0.5), size=3, vjust=-0.6, hjust=-0.5) +
  geom_errorbar(aes(ymin=lower, ymax=upper), size= 0.9, width=0.2, position = position_dodge(width = 0.5)) +
  geom_hline(yintercept = 0, colour="grey", linetype = "longdash", size=0.8)+
  theme(axis.text.x=element_text(face="bold"),
        axis.text.y=element_text(face="bold"),
        axis.title=element_text(size=12,face="bold"),
        strip.text.y = element_text(hjust=0,vjust = 1,angle=180,face="bold")) +
  theme_bw()+
  scale_color_brewer(palette = "Set1")+
  theme(panel.border = element_blank(),
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank(),  
  axis.line.x = element_line(colour = "black"),
  axis.line.y= element_blank()) 

#ggsave(file="lrr.svg", width=9, height=5)


#ggplot(models, aes(paradigm, lrr, fill = intervention)) +
 # ylim(c(-2,2)) +
#  geom_point(shape = 23, size = 2, position = position_dodge(width = 0.5)) + 
#  labs(x = "", y = "lrr", fill = "") +
#  coord_flip() +
#  geom_hline(yintercept = 0, colour="grey", linetype = "longdash")+
#  theme(axis.text.x=element_text(face="bold"),
#  axis.text.y=element_text(face="bold"),
#  axis.title=element_text(size=12,face="bold"),
#  strip.text.y = element_text(hjust=0,vjust = 1,angle=180,face="bold"))+
#  scale_color_brewer(palette = "Set1")

Interpretations

1. Specific techniques are infrequently replicated in drylands globally.
   
2. Intervention and outcome classification useful and appropriate.
   
3. Active versus passive interesting and fascinating.
   
4. Use only random effect model statistics because of heterogeneity and conceptual focus of data and synthesis.
   
5. Active and passive are not the same - active is net positive and passive is net negative. Amazing finding.
6. Then subset out each of active passive and run a meta. Found that both for active and passive, vegetation is a path forward for restoration (list techniques in paper) particularly active.
7. For active interventions aridity index and experiment length are significants in explaining heterogeneity of outcomes. While for passive interventions none of these two moderators were significants.
8. Soil remediation is powerful active intervention and ignoring passive changes in soil is a dramatic and negative impediment to restoration.
   
9. Habitat, animals and soils are challenging to actively (and passively) restore but plants are viable and significant positive outcomes that can be restored.
   
10. Passive restoration confirms that plants can recover passively, but soil does not recover.